############################################################################
# SCRIPT: Municipal-level analysis
############################################################################

################################################################################
##### Packages and directories #####
################################################################################

## Add your own root directory
#setwd("Replication package/Municipal")


## Clean the workspace
rm(list=ls())

## Load packages
library(tidyverse)
library(stargazer)
library(interflex)
library(fixest)
library(ggpubr)
library(fwildclusterboot)
library(xtable)
library(tidylog)
library(sf)
library(tmap)

## Directory of the data
indir <- "Data/"

## Directory of the shapefiles 
inshpdir <- "Data/Shapefiles/"

## Directory of the code
scrdir <- "Code/"

## Directory where to save tables
outtabledir <- "Output/Tables/"

## Directory where to save figures
outplotdir <- "Output/Figures/"



## Functions used in the analysis
source(paste0(scrdir, "functions.R"))


################################################################################
##### Data preparation #####
################################################################################
## Load dataset for analysis
load(paste0(indir, "data_final.RData"))

##### Rescale the signal strength variables #####

## Convert from dBm to dBuV
data_final$strength_loy <- data_final$strength_loy + 107
data_final$freespacestrength_loy <- data_final$freespacestrength_loy + 107

data_final$strength_rpss <- data_final$strength_rpss + 107
data_final$freespacestrength_rpss <- data_final$freespacestrength_rpss + 107


##### Compute the average of signal strength for RPL and RPS #####
data_final$avg_strength_gip <- (data_final$strength_loy + data_final$strength_rpss)/2
data_final$avg_freespacestrength_gip <- (data_final$freespacestrength_loy + data_final$freespacestrength_rpss)/2

##### Dummy for any civilian victim #####
data_final$any_civ_victim <- ifelse(data_final$civ_victims > 0, 1, data_final$civ_victims)


##### Indicator for provincial capital #####
data_final$capital <- ifelse(data_final$mun %in% c("Bilbao", "Vitoria-Gasteiz", "Donostia/San Sebastián", "Pamplona/Iruña"), 1, 0)


##### Compute vote shares over eligible voters (to clean from turnout effects)   #####
# Euskadi: PNV vote share in 1977
data_final$pnv_share_77<- data_final$pnv_77/data_final$censo_77

# Euskadi: Abertzale vote shares in 1977
data_final$ee_share_77 <- data_final$ee_77/data_final$censo_77

data_final$esb_share_77 <- data_final$esb_77/data_final$censo_77

data_final$anv_share_77 <- data_final$anv_77/data_final$censo_77

# --> for Alava, treat Ind14 as EE (following Clark and Llera Ramo, Linz considers them "Independent left")
data_final$ind14_share_77 <- data_final$ind14_77/data_final$censo_77

# Total of Abertzale vote shares in Euskadi
data_final$IA_share_77 <- rowSums(data_final[,c("ee_share_77", "esb_share_77", "anv_share_77", "ind14_share_77")], na.rm = T)

data_final[data_final$mun == "Lantarón",]$IA_share_77 <- NA

## --> Takes value 0 in Navarra


# Navarra: UNAI vote share in 1977
data_final$unai_share_77 <- data_final$unai_77/data_final$censo_77

# Now use UNAI for IA in Navarra
data_final$IA_share_77 <- ifelse(data_final$prov == "Navarra", data_final$unai_share_77, data_final$IA_share_77)

## PSOE vote share in 1977
data_final$psoe_share_77 <- data_final$psoe_77/data_final$censo_77

## NO vote share in 1978 
data_final$no_share_78 <- data_final$no_78/data_final$censo_78

## PNV vote share in 1979
data_final$pnv_share_79 <- data_final$pnv_79/data_final$censo_79

## HB vote share in 1979
data_final$hb_share_79 <- data_final$hb_79/data_final$censo_79

## PSOE vote share in 1979
data_final <- data_final %>% mutate(psoe_share_79 = psoe_79/censo_79)

## UCD vote share in 1979
data_final <- data_final %>% mutate(ucd_share_79 = ucd_79/censo_79)


##### Average vote shares in national elections #####
## PNV
data_final <- data_final %>% 
  mutate(pnv_share_avg = rowMeans(select(data_final, c(pnv_share_77, pnv_share_79,
                                                       pnv_share_82, pnv_share_86,
                                                       pnv_share_89, pnv_share_93,
                                                       pnv_share_96, pnv_share_00,
                                                       pnv_share_04, pnv_share_08,
                                                       pnv_share_11, pnv_share_15, 
                                                       pnv_share_16, pnv_share_19)),
                                  na.rm = T))

## Radical nationalists (Izquierda Abertzale): IA/HB
data_final <- data_final %>% 
  mutate(ia_share_avg = rowMeans(select(data_final, c(IA_share_77, hb_share_79,
                                                      hb_share_82, hb_share_86,
                                                      hb_share_89, hb_share_93,
                                                      hb_share_96, amaiur_share_11,
                                                      bildu_share_15, bildu_share_16,
                                                      bildu_share_19)),
                                 na.rm = T))

## PSOE
data_final <- data_final %>%
  mutate(psoe_share_avg = rowMeans(select(data_final, c(psoe_share_77, psoe_share_79,
                                                        psoe_share_82, psoe_share_86,
                                                        psoe_share_89, psoe_share_93,
                                                        psoe_share_96, psoe_share_00,
                                                        psoe_share_04, psoe_share_08,
                                                        psoe_share_11, psoe_share_15,
                                                        psoe_share_16, psoe_share_19)),
                                   na.rm = T))

## CP/PP
data_final <- data_final %>%
  mutate(pp_share_avg = rowMeans(select(data_final, c(cp_share_82, cp_share_86,
                                                      pp_share_89, pp_share_93,
                                                      pp_share_96, pp_share_00,
                                                      pp_share_04, pp_share_08,
                                                      pp_share_11, pp_share_15,
                                                      pp_share_16, pp_share_19)),
                                 na.rm = T))


##### Average vote shares in regional elections #####
## Note that the BC and Navarra have different election years

## PNV (only for BC)
data_final <- data_final %>%
  mutate(pnv_er_share_avg=rowMeans(select(data_final, c(pnv_share_er80, pnv_share_er84,
                                                        pnv_share_er86, pnv_share_er90,
                                                        pnv_share_er94, pnv_share_er98,
                                                        pnv_share_er01, pnv_share_er05,
                                                        pnv_share_er09, pnv_share_er12,
                                                        pnv_share_er16, pnv_share_er20)),
                                   na.rm=T))


## IA (both BC and Navarre)
data_final <- data_final %>%
  mutate(ia_er_share_avg=rowMeans(select(data_final, c(hb_share_er79, hb_share_er80,
                                                       hb_share_er83, hb_share_er84, 
                                                       hb_share_er86, hb_share_er87,
                                                       hb_share_er90, hb_share_er91,
                                                       hb_share_er94, hb_share_er95,
                                                       eh_share_er98, eh_share_er99,
                                                       eh_share_er01, bildu_share_er11,
                                                       bildu_share_er12, bildu_share_er15,
                                                       bildu_share_er16, bildu_share_er19,
                                                       bildu_share_er20, bildu_share_er23)),
                                  na.rm=T))



## PSOE (only BC)
data_final <- data_final %>%
  mutate(psoe_er_share_avg=rowMeans(select(data_final, c(psoe_share_er80, psoe_share_er84,
                                                         psoe_share_er86, psoe_share_er90,
                                                         psoe_share_er94, psoe_share_er98,
                                                         psoe_share_er01, psoe_share_er05,
                                                         psoe_share_er09, psoe_share_er12,
                                                         psoe_share_er16, psoe_share_er20)),
                                    na.rm=T))

## AP/CP/PP (only BC)
data_final <- data_final %>%
  mutate(pp_er_share_avg=rowMeans(select(data_final, c(ap_share_er80, ap_share_er84,
                                                       ap_share_er86, pp_share_er90,
                                                       pp_share_er94, pp_share_er98,
                                                       pp_share_er01, pp_share_er05,
                                                       pp_share_er09, pp_share_er12,
                                                       pp_share_er16, pp_share_er20)),
                                  na.rm=T))



## Give missing value to towns dissolved in other towns
#  after 1977 (to avoid double-counting them)
dissolved_post_77 <- c("Arlucea-Marquínez", "Barriobusto", "Bergüenda", 
                       "Labraza", "Pipaón", "Salcedo")

data_final$pnv_share_avg[data_final$mun%in%dissolved_post_77] <- NA
data_final$ia_share_avg[data_final$mun%in%dissolved_post_77] <- NA
data_final$pp_share_avg[data_final$mun%in%dissolved_post_77] <- NA
data_final$psoe_share_avg[data_final$mun%in%dissolved_post_77] <- NA

data_final$pnv_er_share_avg[data_final$mun%in%dissolved_post_77] <- NA
data_final$ia_er_share_avg[data_final$mun%in%dissolved_post_77] <- NA
data_final$psoe_er_share_avg[data_final$mun%in%dissolved_post_77] <- NA
data_final$pp_er_share_avg[data_final$mun%in%dissolved_post_77] <- NA

## Compute changes in basque speakers over time
data_final <- data_final %>% mutate(change_72_16 = eusk_share_16 - euskaldun_share)



##### Compute victims of violence per capita #####
data_final$civ_victims_pc <- data_final$civ_victims/data_final$pop_30
data_final$civ_victims_pc[data_final$mun == "Lezáun"] <- 0

data_final$shot_pc <- data_final$shot/data_final$pop_30
data_final$shot_pc[data_final$mun == "Lezáun"] <- 0


################################################################################
##### Figure 2 and Figure 3: Maps #####
################################################################################
## Import the shapefile of 1970 municipal borders
bd <- read_sf(dsn=paste0(inshpdir, "boundaries.shp"))

## Rename variable names to match data 
bd <- rename(bd, mun = NMUN, prov = NPRO)

## Harmonization of provincial and municipal names
bd <- bd %>% mutate(prov = str_replace_all(prov, "Araba/Álava", "Araba-Alava")) %>%
  mutate(mun = str_replace_all(mun, " de ", " De ")) %>%
  mutate(mun = str_replace_all(mun, " el ", " El "))


bd$mun[bd$mun=="Arrazua-Ubarrundia"] <- "Arratzua-Ubarrundia"
bd$mun[bd$mun=="Arrazua-Ubarrundia"] <- "Arratzua-Ubarrundia"
bd$mun[bd$mun=="Moreda De Álava"] <- "Moreda De Álava/Moreda Araba"
bd$mun[bd$mun=="Salvatierra/Agurain"] <- "Agurain/Salvatierra"
bd$mun[bd$mun=="Bidegoian"] <- "Bidania-Goiatz"
bd$mun[bd$mun=="Arrasate/Mondragón"] <- "Arrasate/Mondragon"
bd$mun[bd$mun=="Soraluze/Placencia De las Armas"] <- "Soraluze-Placencia De Las Armas"
bd$mun[bd$mun=="Donostia-San Sebastián"] <- "Donostia/San Sebastián"
bd$mun[bd$mun=="Abárzuza"] <- "Abárzuza/Abartzuza"
bd$mun[bd$mun=="Allín"] <- "Allín/Allin"
bd$mun[bd$mun=="Arcos, Los"] <- "Los Arcos"
bd$mun[bd$mun=="Atez"] <- "Atez/Atetz"
bd$mun[bd$mun=="Busto, El"] <- "El Busto"
bd$mun[bd$mun=="Ciriza"] <- "Ciriza/Ziritza"
bd$mun[bd$mun=="Etxarri-Aranatz"] <- "Etxarri Aranatz"
bd$mun[bd$mun=="Estella/Lizarra"] <- "Estella-Lizarra"
bd$mun[bd$mun=="Leache"] <- "Leache/Leatxe"
bd$mun[bd$mun=="Monreal"] <- "Monreal/Elo"
bd$mun[bd$mun=="Navascués"] <- "Navascués/Nabaskoze"
bd$mun[bd$mun=="Orbaitzeta"] <- "Orbaizeta"
bd$mun[bd$mun=="Puente la Reina/Gares"] <- "Puente La Reina/Gares"
bd$mun[bd$mun=="Saldías"] <- "Saldias"
bd$mun[bd$mun=="Torralba del Río"] <- "Torralba Del Río"
bd$mun[bd$mun=="Torres del Río"] <- "Torres Del Río"
bd$mun[bd$mun=="Ucar"] <- "Úcar"
bd$mun[bd$mun=="Unzué"] <- "Unzué/Untzue"
bd$mun[bd$mun=="Urraul Alto"] <- "Urraúl Alto"
bd$mun[bd$mun=="Urraul Bajo"] <- "Urraúl Bajo"
bd$mun[bd$mun=="Abanto y Ciérvana-Abanto Zierbena"] <- "Abanto Y Ciérvana-Abanto Zierbena"
bd$mun[bd$mun=="Sopelana"] <- "Sopela"

## See municipalities not present in the map
anti_join(bd, data_final, by=c("prov", "mun")) %>% st_drop_geometry() %>% select(c(prov, mun)) %>% print(n=Inf)
anti_join(data_final, bd, by=c("prov", "mun")) %>% select(c(prov, mun))

## Merge variables to the polygons
bd <- merge(bd, data_final, by=c("prov", "mun"), all.x=T)

## Coordinates of the RPL transmitter (for plotting)
antenna <- data_final[1, c("lon_loy", "lat_loy")]

## Convert to the CRS of the polygons (to draw the map)
antenna <- st_as_sf(antenna, coords=c("lon_loy", "lat_loy"), crs="+proj=longlat")
antenna <- st_transform(antenna, crs=st_crs(bd))


## Map of RPL signal strength
tm_shape(filter(bd, eus==1)) +
  tm_borders() +
  tm_fill(col="strength_loy", n=10, title="Signal strength", midpoint=NA, palette = "viridis") +
  tm_layout(legend.outside = T, 
            legend.outside.position = "right",
            legend.outside.size = 0.3) +
  tm_shape(antenna) + 
  tm_bubbles(size=0.1, shape = 17, col = "black") +
  tm_add_legend("symbol", col = "black", shape=17, size = c(.5), labels = "Transmitter")


## Save and export
tmap_save(filename = paste0(outplotdir, "figure2.pdf"),
          width = 10, height = 6, units = "cm")


## Map of PNV vote share 1933
tm_shape(filter(bd, eus==1)) +
  tm_borders() +
  tm_fill(col="pnv_share_33", n=10, title="", midpoint=NA, palette = "Greens") +
  tm_layout(legend.outside = T, 
            legend.outside.position = "left",
            legend.outside.size = 0.3,
            legend.text.size = 2)
## Save and export
tmap_save(filename = paste0(outplotdir, "figure3a.pdf"))

## Map of euskera speakers 1972
tm_shape(filter(bd, eus==1)) +
  tm_borders() +
  tm_fill(col="euskaldun_share", n=10, title="", midpoint=NA, palette = "Blues") +
  tm_layout(legend.outside = T, 
            legend.outside.position = "right",
            legend.outside.size = 0.3,
            legend.text.size = 2)
## Save and export
tmap_save(filename = paste0(outplotdir, "figure3b.pdf"))

## Clean the memory
rm(antenna)

################################################################################
##### Further data preparation #####
################################################################################


## Normalize signal strength strength to be [0,1]
data_final$strength_loy <- (data_final$strength_loy - min(data_final$strength_loy))/(max(data_final$strength_loy) - min(data_final$strength_loy))
data_final$strength_rpss <- (data_final$strength_rpss - min(data_final$strength_rpss))/(max(data_final$strength_rpss) - min(data_final$strength_rpss))
data_final$avg_strength_gip <- (data_final$avg_strength_gip - min(data_final$avg_strength_gip))/(max(data_final$avg_strength_gip) - min(data_final$avg_strength_gip))
  

##### Create additional variables for analysis #####
## Log pop 1970
data_final$log_pop_70 <- log(data_final$pop_70)

# Log pop 1960 (for balance only)
data_final$log_pop_60 <- log(data_final$pop_60)

## Log area
data_final$log_area <- log(data_final$area)

## Log height
data_final$log_height <- log(data_final$height)

## Log ETA prisoners
data_final$log_pris <- log(1 + data_final$eta_pris)

## Any ETA prisoner
data_final$any_pris <- ifelse(data_final$eta_pris>0,1,0)

## ETA attacks and victims (binary indicators)
data_final <- data_final %>% mutate(any_eta_att_68_77 = case_when(eta_att_68_77 > 0 ~ 1,
                                                                  eta_att_68_77 == 0 ~ 0),
                                    any_eta_att_78_94 = case_when(eta_att_78_94 > 0 ~ 1,
                                                                  eta_att_78_94 == 0 ~ 0),
                                    any_eta_att_95_03 = case_when(eta_att_95_03 > 0 ~ 1,
                                                                  eta_att_95_03 == 0 ~ 0),
                                    any_eta_vict_68_77 = case_when(eta_vict_68_77 > 0 ~ 1,
                                                                   eta_vict_68_77 == 0 ~ 0),
                                    any_eta_vict_78_94 = case_when(eta_vict_78_94 > 0 ~ 1,
                                                                   eta_vict_78_94 == 0 ~ 0),
                                    any_eta_vict_95_03 = case_when(eta_vict_95_03 > 0 ~ 1,
                                                                   eta_vict_95_03 == 0 ~ 0))

## Log ETA attacks and victims
data_final <- data_final %>% mutate(log_eta_att_68_77 = log(1+eta_att_68_77),
                                    log_eta_att_78_94 = log(1+eta_att_78_94),
                                    log_eta_att_95_03 = log(1+eta_att_95_03),
                                    log_eta_vict_68_77 = log(1+eta_vict_68_77),
                                    log_eta_vict_78_94 = log(1+eta_vict_78_94),
                                    log_eta_vict_95_03 = log(1+eta_vict_95_03))

## Distance from provincial capitals
data_final$dist_capital <- case_when(data_final$prov == "Araba-Alava" ~ data_final$dist_vi,
                                     data_final$prov == "Bizkaia" ~ data_final$dist_bi,
                                     data_final$prov == "Gipuzkoa" ~ data_final$dist_ss,
                                     data_final$prov == "Navarra" ~ data_final$dist_pamp)

## Log distance from provincial capitals
data_final$log_dist_capital <- log(1 + data_final$dist_capital)

## Rename the FS variable
data_final <- data_final %>% rename(fs_loy = freespacestrength_loy)


################################################################################
##### Model and table settings #####
################################################################################

## Set Conley covariance matrix (30 km cutoff)
conley_vcov <- fixest::vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 30)

## Dictionary to label variables in tables
setFixest_dict(c(strength_loy = "Signal strength", pnv_share_avg = "PNV (Avg)",
                 ia_share_avg = "Rad. Nat. (Avg)", pp_share_avg = "PP (Avg)",
                 psoe_share_avg = "PSOE (Avg)", 
                 pnv_er_share_avg = "PNV (Reg. Avg.)", ia_er_share_avg = "Rad. Nat. (Reg. Avg.)",
                 pp_er_share_avg = "PP (Reg. Avg.)", psoe_er_share_avg = "PSOE (Reg. Avg.)",
                 "strength_loy:pnv_share_33" = "Signal strength * PNV '33",
                 "strength_loy:euskaldun_share" = "Signal strength * Sh. Basque speakers",
                 pnv_share_33 = "PNV '33", euskaldun_share = "Sh. Basque speakers",
                 log_pris = "ETA pris.", turnout_78 = "Turnout '78", 
                 no_share_78 = "No share '78", turnout_r79 = "Turnout '79",
                 yes_share_r79 = "Yes share '79", change_72_16 = "Delta speakers (1972-2016)",
                 log_eta_att_68_77 = "Log(1+Att.) 68-77", log_eta_att_78_94 = "Log(1+Att.) 78-94", 
                 log_eta_att_95_03 = "Log(1+Att.) 95-03", any_eta_att_68_77="Att.>0 68-77", 
                 any_eta_att_78_94="Att.>0 78-94", any_eta_att_95_03="Att.>0 95-03", 
                 log_eta_vict_68_77 = "Log(1+Vict.) 68-77", log_eta_vict_78_94 = "Log(1+Vict.) 78-94", log_eta_vict_95_03 = "Log(1+Vict.) 95-03",
                 any_eta_vict_68_77 = "Vict.>0 68-77", any_eta_vict_78_94 = "Vict.>0, 78-94", any_eta_vict_95_03 = "Vict.>0, 95-03",
                 avg_strength_gip="Avg. RPL + RPS", strength_rpss="Signal strength (RPS)", ucd_share_79="UCD (1979)"))



################################################################################
##### Table A1: Summary statistics #####
################################################################################
data_final %>% filter(eus==1) %>% select(c(strength_loy, fs_loy, log_height,
                                           log_area, log_pop_70, log_pop_60, ia_share_avg,
                                           pnv_share_avg, pp_share_avg, 
                                           psoe_share_avg, turnout_78, turnout_r79,
                                           no_share_78, yes_share_r79,
                                           eta_pris, any_civ_victim, civ_victims_pc,
                                           shot_pc, pnv_share_33, left_share_33,
                                           right_share_33, yes_share_33,
                                           euskaldun_share, log_dist_capital,
                                           change_72_16, jm_33, batzoki_33,
                                           mend_33, eab_36, 
                                           eta_att_68_77, eta_att_78_94, eta_att_95_03,
                                           eta_vict_68_77, eta_vict_78_94, eta_vict_95_03
                                           )) %>%
  stargazer(out = paste0(outtabledir, "tableA1.tex"),
            covariate.labels = c("Signal strength (stand.)", "Signal strength in f.s.", 
                                 "Log(height)", "Log(area)", 
                                 "Log(Pop. 1970)", "Log(Pop. 1960)", "Rad. nat. vote share (Avg)",
                                 "PNV vote share (Avg)", "PP vote share (Avg)", 
                                 "PSOE vote share (Avg)", "Turnout ref. 1978", 
                                 "Turnout ref. 1979", "No share ref. 1978", 
                                 "Yes share ref. 1979", "Eta prisoners (raw)", 
                                 "Civilian victims (dummy)", "Civilian victims p.c.",
                                 "Executions p.c.", "PNV vote share 1933", 
                                 "Left vote share 1933", "Right vote share 1933",
                                 "Yes share ref. 1933", "Basque speakers share 1972",
                                 "Log(dist. capital)", "Delta Basque speakers 1972-2016",
                                 "PNV branch 1933", "Batzoki 1933", "Alpinist group 1933",
                                 "Patriotic Women 1936", "ETA attacks 68-77 (raw)", "ETA attacks 78-94 (raw)",
                                 "ETA attacks 95-03 (raw)", "ETA victims 68-77 (raw)",
                                 "ETA victims 78-94 (raw)", "ETA victims 95-03 (raw)"
                                 ))
 
################################################################################
##### Table 1: Balance tests #####
################################################################################

########## Balance tests for Euskadi ##########

## Names of variables used in the balance tests
bnames <- c("Civ. victims  [0/1]", "Civ. victims p.c.", "Executions p.c.", 
            "% PNV  1933", "% Left  1933", "% Right  1933", "% Yes  1933", 
            "Basque  speakers 1972", "Log(Dist. capital)", "Log(Pop. 1960)",  "PNV local branch 1933 [0/1]",
            "Batzoki 1933 [0/1]", "Basque Youth 1933 [0/1]", "Alpinism group 1933 [0/1]",
            "Patriotic Women 1936 [0/1]")

##### Without PNV '33 #####
btests <- NULL

# Any civilian victim
btests[[1]] <- feols(any_civ_victim ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)

# Civilian victims p.c.
btests[[2]] <- update(btests[[1]], scale(civ_victims_pc) ~ .)

# Executions p.c.
btests[[3]] <- update(btests[[1]], scale(shot_pc) ~ .)

# % PNV 1933
btests[[4]] <- update(btests[[1]], scale(pnv_share_33) ~ .)

# % Left 1933
btests[[5]] <- update(btests[[1]], scale(left_share_33) ~ .)

# % Right 1933
btests[[6]] <- update(btests[[1]], scale(right_share_33) ~ .)

# % Yes 1933
btests[[7]] <- update(btests[[1]], scale(yes_share_33) ~ .)

# % Basque speakers 1972
btests[[8]] <- update(btests[[1]], scale(euskaldun_share) ~ .)

# Log distance from prov. capital
btests[[9]] <- update(btests[[1]], scale(log_dist_capital) ~ .)

# Log pop. 1960
btests[[10]] <- update(btests[[1]], scale(log_pop_60) ~ .)

# PNV juntas municipales 1933
btests[[11]] <- update(btests[[1]], jm_33 ~ .)

## Batzoki 1933
btests[[12]] <- update(btests[[1]], batzoki_33 ~ .)

## Juventud Vasca 1933
btests[[13]] <- update(btests[[1]], jv_33 ~ .)

## Mendigotzale 1933
btests[[14]] <- update(btests[[1]], mend_33 ~ .)

## EAB 1936
btests[[15]] <- update(btests[[1]], eab_36 ~ .)


##### Controlling for PNV '33 #####
btests_augm <- NULL

# Any civilian victim
btests_augm[[1]] <- update(btests[[1]], . ~ . + pnv_share_33)

# Civilian victims p.c.
btests_augm[[2]] <- update(btests_augm[[1]], scale(civ_victims_pc) ~ .)

# Executions p.c.
btests_augm[[3]] <- update(btests_augm[[1]], scale(shot_pc) ~ .)

# % Left 1933
btests_augm[[4]] <- update(btests_augm[[1]], scale(left_share_33) ~ .)

# % Right 1933
btests_augm[[5]] <- update(btests_augm[[1]], scale(right_share_33) ~ .)

# % Yes 1933
btests_augm[[6]] <- update(btests_augm[[1]], scale(yes_share_33) ~ .)

# % Basque speakers 1972
btests_augm[[7]] <- update(btests_augm[[1]], scale(euskaldun_share) ~ .)

# Log distance from prov. capital
btests_augm[[8]] <- update(btests_augm[[1]], scale(log_dist_capital) ~  .)

# Log pop. 1960
btests_augm[[9]] <- update(btests_augm[[1]], scale(log_pop_60) ~ .)

# PNV juntas municipales 1933
btests_augm[[10]] <- update(btests_augm[[1]], jm_33 ~ .)

## Batzoki 1933
btests_augm[[11]] <- update(btests_augm[[1]], batzoki_33 ~ .)

## Juventud Vasca 1933
btests_augm[[12]] <- update(btests_augm[[1]], jv_33 ~ .)

## Mendigotzale 1933
btests_augm[[13]] <- update(btests_augm[[1]], mend_33 ~ .)

## EAB 1936
btests_augm[[14]] <- update(btests_augm[[1]], eab_36 ~ .)


## Extract statistics for balance table, customize, and print
b <- do.call("rbind", lapply(btests, balance_coefs))

bb <- rbind(do.call("rbind", lapply(btests_augm[1:3], balance_coefs)),
            cbind("pnv_share_33", t(rep("-", 3))),
            do.call("rbind", lapply(btests_augm[4:14], balance_coefs)))

btable <- as.matrix(cbind(b, bb, do.call("rbind", lapply(btests, nobs))))[,-5]
btable[,"depvar"] <- bnames
colnames(btable) <- c("Dep. Var.", "Coef.", "SE", "p", "Coef.", "SE", "p", "N")

print(xtable(btable, caption = "Balance tests", label = "tab:balance_table", align = "llccccccc"),
      add.to.row = list(pos = list(0,0,0),
                        command = c("& \\multicolumn{3}{c}{Without PNV} & \\multicolumn{3}{c}{Including PNV} & \\\\",
                                    "\\cmidrule(lr){2-4} \\cmidrule(lr){5-7} \\addlinespace",
                                    "Dep. Var. & Coef. & SE & p & Coef. & SE & p & N \\\\")),
      include.colnames = F,
      include.rownames = F,
      caption.placement = "top",
      type = "latex",
      file = paste0(outtabledir, "table1.tex"), 
      replace = TRUE
)


###################################################
##### Descriptive results after balance tests #####
###################################################

## SD of the residualized treatment (mentioned in text)
lm(strength_loy~ fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 + factor(prov), 
   data=subset(data_final, eus==1)) %>% residuals %>% sd()

##### Figure B1 #####

## Correlation between signal and contributors to RPL
l1 <- data_final %>%
  filter(eus==1) %>%
  mutate(subscriptions_pc = subscriptions_78/pop_70) %>%
  ggplot(aes(x=strength_loy, y=subscriptions_pc)) +
  geom_point() +
  labs(x="", y="RPL subscriptions p.c.") + 
  geom_smooth(method="loess") +
  theme_bw()

## Correlation between signal and RPL delegations
l2 <- data_final %>%
  filter(eus==1) %>%
  mutate(rpl_delegation = case_when(delegations_78 > 0 ~ 1,
                                    delegations_78 == 0 ~ 0)) %>%
  ggplot(aes(x=strength_loy, y=rpl_delegation)) +
  geom_point() +
  labs(x="", y="RPL delegation") +
  geom_smooth(method="loess") +
  theme_bw()

## Put in unique plot and export
ggarrange(l1, l2, nrow = 1) %>% annotate_figure(bottom = "RPL signal strength")

ggsave(filename = "figureB1.pdf", 
       path = outplotdir,
       width = 12, units = "cm")


################################################################################
##### Table 2: Average effects on voting #####
################################################################################


##### Average vote shares in the post-Franco era #####

### Outcomes ###
## PNV: 1977-2019
## IA: IA total (1977), HB (1979-1996), Amaiur (2011), EH Bildu (2015-2019)
## PP: CP (1982-1986), PP (1989-2019)
## PSOE: 1977-2019

## Regressions on average outcomes
lt_avgs <- NULL

## IA
lt_avgs[[1]] <- feols(ia_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                      data = subset(data_final, eus==1),
                      vcov = conley_vcov)

## PNV
lt_avgs[[2]] <- update(lt_avgs[[1]], pnv_share_avg ~ .)

## PP
lt_avgs[[3]] <- update(lt_avgs[[1]], pp_share_avg ~ .)

## PSOE
lt_avgs[[4]] <- update(lt_avgs[[1]], psoe_share_avg ~ .)

## Export table
etable(lt_avgs, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "table2.tex"), replace = TRUE)


################################################################################
##### Figure 4: Results in single elections #####
################################################################################

# Post-Franco general election years
election_years <- c("1977", "1979", "1982", "1986", "1989", "1993", "1996",
                    "2000", "2004", "2008", "2011", "2015", "2016", "2019")


## PNV
pnv_yby <- NULL

pnv_yby[[1]] <- feols(pnv_share_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                      data = subset(data_final, eus==1),
                      vcov = conley_vcov)
pnv_yby[[2]] <- update(pnv_yby[[1]], pnv_share_79 ~ .)
pnv_yby[[3]] <- update(pnv_yby[[1]], pnv_share_82 ~ .)
pnv_yby[[4]] <- update(pnv_yby[[1]], pnv_share_86 ~ .)
pnv_yby[[5]] <- update(pnv_yby[[1]], pnv_share_89 ~ .)
pnv_yby[[6]] <- update(pnv_yby[[1]], pnv_share_93 ~ .)
pnv_yby[[7]] <- update(pnv_yby[[1]], pnv_share_96 ~ .)
pnv_yby[[8]] <- update(pnv_yby[[1]], pnv_share_00 ~ .)
pnv_yby[[9]] <- update(pnv_yby[[1]], pnv_share_04 ~ .)
pnv_yby[[10]] <- update(pnv_yby[[1]], pnv_share_08 ~ .)
pnv_yby[[11]] <- update(pnv_yby[[1]], pnv_share_11 ~ .)
pnv_yby[[12]] <- update(pnv_yby[[1]], pnv_share_15 ~ .)
pnv_yby[[13]] <- update(pnv_yby[[1]], pnv_share_16 ~ .)
pnv_yby[[14]] <- update(pnv_yby[[1]], pnv_share_19 ~ .)


## Store coefficients and CIs in a matrix and plot
pnv_plot <- 
  cbind(
    # IDs
    c(1:length(pnv_yby)),
    # Years
    as.factor(election_years),
    # Coefficients
    do.call("rbind", lapply(pnv_yby, function(x) coef(x)["strength_loy"])),
    # 95% CIs
    do.call("rbind", lapply(pnv_yby, function(x) confint(x)["strength_loy",]))
  ) %>%
  # Turn into df
  as.data.frame() %>%
  dplyr::rename(year = `c(1:length(pnv_yby))`, coefs =  strength_loy, lower = `2.5 %`, upper = `97.5 %`) %>%
  # Plot
  ggplot(aes(x=year, y = coefs)) + geom_point() + geom_errorbar(aes(ymin = lower, ymax = upper), color="black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "blue") +
  scale_x_discrete(limits = 1:length(pnv_yby), labels = election_years, drop=F) +
  labs(x = "", y = "", title = "PNV") + theme_bw() + theme(axis.text.x = element_text(size=6.5, angle = 45))


## IA
ia_yby <- NULL

ia_yby[[1]] <- feols(IA_share_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)
ia_yby[[2]] <- update(ia_yby[[1]], hb_share_79 ~ .)
ia_yby[[3]] <- update(ia_yby[[1]], hb_share_82 ~ .)
ia_yby[[4]] <- update(ia_yby[[1]], hb_share_86 ~ .)
ia_yby[[5]] <- update(ia_yby[[1]], hb_share_89 ~ .)
ia_yby[[6]] <- update(ia_yby[[1]], hb_share_93 ~ .)
ia_yby[[7]] <- update(ia_yby[[1]], hb_share_96 ~ .)
ia_yby[[8]] <- update(ia_yby[[1]], amaiur_share_11 ~ .)
ia_yby[[9]] <- update(ia_yby[[1]], bildu_share_15 ~ .)
ia_yby[[10]] <- update(ia_yby[[1]], bildu_share_16 ~ .)
ia_yby[[11]] <- update(ia_yby[[1]], bildu_share_19 ~ .)

ia_plot <- 
  cbind(
    # IDs
    c(1:length(ia_yby)),
    # Years
    as.factor(election_years[!election_years%in%c("2000", "2004", "2008")]),
    # Coefficients
    do.call("rbind", lapply(ia_yby, function(x) coef(x)["strength_loy"])),
    # 95% CIs
    do.call("rbind", lapply(ia_yby, function(x) confint(x)["strength_loy",]))
  ) %>%
  # Turn into df
  as.data.frame() %>%
  dplyr::rename(year = `c(1:length(ia_yby))`, coefs =  strength_loy, lower = `2.5 %`, upper = `97.5 %`) %>%
  # Plot
  ggplot(aes(x=year, y = coefs)) + geom_point() + geom_errorbar(aes(ymin = lower, ymax = upper), color="black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "blue") +
  scale_x_discrete(limits = 1:length(ia_yby), labels = election_years[!election_years%in%c("2000", "2004", "2008")], drop=F) +
  labs(x = "", y = "", title = "Radical nationalists") + theme_bw() + theme(axis.text.x = element_text(size=6.5, angle = 45))



## PP
pp_yby <- NULL

pp_yby[[1]] <- feols(cp_share_82 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)
pp_yby[[2]] <- update(pp_yby[[1]], cp_share_86 ~ .)
pp_yby[[3]] <- update(pp_yby[[1]], pp_share_89 ~ .)
pp_yby[[4]] <- update(pp_yby[[1]], pp_share_93 ~ .)
pp_yby[[5]] <- update(pp_yby[[1]], pp_share_96 ~ .)
pp_yby[[6]] <- update(pp_yby[[1]], pp_share_00 ~ .)
pp_yby[[7]] <- update(pp_yby[[1]], pp_share_04 ~ .)
pp_yby[[8]] <- update(pp_yby[[1]], pp_share_08 ~ .)
pp_yby[[9]] <- update(pp_yby[[1]], pp_share_11 ~ .)
pp_yby[[10]] <- update(pp_yby[[1]], pp_share_15 ~ .)
pp_yby[[11]] <- update(pp_yby[[1]], pp_share_16 ~ .)
pp_yby[[12]] <- update(pp_yby[[1]], pp_share_19 ~ .)

pp_plot <- 
  cbind(
    # IDs
    c(1:length(pp_yby)),
    # Years
    as.factor(election_years[!election_years%in%c("1977", "1979")]),
    # Coefficients
    do.call("rbind", lapply(pp_yby, function(x) coef(x)["strength_loy"])),
    # 95% CIs
    do.call("rbind", lapply(pp_yby, function(x) confint(x)["strength_loy",]))
  ) %>%
  # Turn into df
  as.data.frame() %>%
  dplyr::rename(year = `c(1:length(pp_yby))`, coefs =  strength_loy, lower = `2.5 %`, upper = `97.5 %`) %>%
  # Plot
  ggplot(aes(x=year, y = coefs)) + geom_point() + geom_errorbar(aes(ymin = lower, ymax = upper), color="black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "blue") +
  scale_x_discrete(limits = 1:length(pp_yby), labels = election_years[!election_years%in%c("1977", "1979")], drop=F) +
  labs(x = "", y = "", title = "PP") + theme_bw() + theme(axis.text.x = element_text(size=6.5, angle = 45))


## PSOE: 1977-2019
psoe_yby <- NULL

psoe_yby[[1]] <- feols(psoe_share_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                       data = subset(data_final, eus==1),
                       vcov = conley_vcov)
psoe_yby[[2]] <- update(psoe_yby[[1]], psoe_share_79 ~ .)
psoe_yby[[3]] <- update(psoe_yby[[1]], psoe_share_82 ~ .)
psoe_yby[[4]] <- update(psoe_yby[[1]], psoe_share_86 ~ .)
psoe_yby[[5]] <- update(psoe_yby[[1]], psoe_share_89 ~ .)
psoe_yby[[6]] <- update(psoe_yby[[1]], psoe_share_93 ~ .)
psoe_yby[[7]] <- update(psoe_yby[[1]], psoe_share_96 ~ .)
psoe_yby[[8]] <- update(psoe_yby[[1]], psoe_share_00 ~ .)
psoe_yby[[9]] <- update(psoe_yby[[1]], psoe_share_04 ~ .)
psoe_yby[[10]] <- update(psoe_yby[[1]], psoe_share_08 ~ .)
psoe_yby[[11]] <- update(psoe_yby[[1]], psoe_share_11 ~ .)
psoe_yby[[12]] <- update(psoe_yby[[1]], psoe_share_15 ~ .)
psoe_yby[[13]] <- update(psoe_yby[[1]], psoe_share_16 ~ .)
psoe_yby[[14]] <- update(psoe_yby[[1]], psoe_share_19 ~ .)


psoe_plot <- 
  cbind(
    # IDs
    c(1:length(psoe_yby)),
    # Years
    as.factor(election_years),
    # Coefficients
    do.call("rbind", lapply(psoe_yby, function(x) coef(x)["strength_loy"])),
    # 95% CIs
    do.call("rbind", lapply(psoe_yby, function(x) confint(x)["strength_loy",]))
  ) %>%
  # Turn into df
  as.data.frame() %>%
  dplyr::rename(year = `c(1:length(psoe_yby))`, coefs =  strength_loy, lower = `2.5 %`, upper = `97.5 %`) %>%
  # Plot
  ggplot(aes(x=year, y = coefs)) + geom_point() + geom_errorbar(aes(ymin = lower, ymax = upper), color="black") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "blue") +
  scale_x_discrete(limits = 1:length(psoe_yby), labels = election_years, drop=F) +
  labs(x = "", y = "", title = "PSOE") + theme_bw() + theme(axis.text.x = element_text(size=6.5, angle = 45))


## Make a unique plot
ggarrange(ia_plot, pnv_plot, pp_plot, psoe_plot, 
          ncol = 2, nrow = 2)

## Export the plot
ggsave(filename = "figure4.pdf", path = outplotdir,width = 14,  height = 14, units = "cm")


################################################################################
##### Robustness checks #####
################################################################################

##### Table C1: Control also for euskaldun share #####
models <- NULL
models[[1]] <- feols(ia_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 + euskaldun_share | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)
models[[2]] <- update(models[[1]], pnv_share_avg ~ .)
models[[3]] <- update(models[[1]], pp_share_avg ~ .)
models[[4]] <- update(models[[1]], psoe_share_avg ~ .)

## Export table
etable(models, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC1.tex"), replace = TRUE)


##### Table C2: Control for all variables in the balance checks #####
models <- NULL
models[[1]] <- feols(ia_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 +
        civ_victims_pc + shot_pc + left_share_33 + right_share_33 + yes_share_33 +
        euskaldun_share + log_dist_capital + log_pop_60 + jm_33 + batzoki_33 + jv_33 + mend_33 + eab_36 | prov,
        data = subset(data_final, eus==1),
        vcov = conley_vcov)
models[[2]] <- update(models[[1]], pnv_share_avg ~ .)
models[[3]] <- update(models[[1]], pp_share_avg ~ .)
models[[4]] <- update(models[[1]], psoe_share_avg ~ .)

## Export table
etable(models, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC2.tex"), replace = TRUE)



##### Table C3: Replace free space signal with log(lat) and log(lon) #####
models <- NULL
models[[1]] <- feols(ia_share_avg ~ strength_loy + log(abs(lat_rec)) + log(abs(lon_rec)) + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)
models[[2]] <- update(models[[1]], pnv_share_avg ~ .)
models[[3]] <- update(models[[1]], pp_share_avg ~ .)
models[[4]] <- update(models[[1]], psoe_share_avg ~ .)

## Export table
etable(models, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC3.tex"), replace = TRUE)

##### Table C4: Replace free space signal with log + lat interaction #####
models <- NULL
models[[1]] <- feols(ia_share_avg ~ strength_loy + lat_rec*lon_rec + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)
models[[2]] <- update(models[[1]], pnv_share_avg ~ .)
models[[3]] <- update(models[[1]], pp_share_avg ~ .)
models[[4]] <- update(models[[1]], psoe_share_avg ~ .)

## Export table
etable(models, keep = "Signal strength", tex=T, style.tex=style.tex("aer"), file=paste0(outtabledir, "tableC4.tex"), replace = TRUE)



##### Table C5: Alternative Conley cutoffs #####

## Matrix to store the results
conleys <- matrix(nrow=4, ncol=4)

## IA
fit <- feols(ia_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
             data = subset(data_final, eus==1), vcov=conley_vcov)
conleys[1,1] <- coef(fit)["strength_loy"]
# 30 km (base)
conleys[1,2] <- se(fit)["strength_loy"]
# 40 km
conleys[1,3] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 40))["strength_loy"]
# 50 km
conleys[1,4] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 50))["strength_loy"]

## PNV
fit <- feols(pnv_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
             data = subset(data_final, eus==1), vcov=conley_vcov)
conleys[2,1] <- coef(fit)["strength_loy"]
# 30 km (base)
conleys[2,2] <- se(fit)["strength_loy"]
# 40 km
conleys[2,3] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 40))["strength_loy"]
# 50 km
conleys[2,4] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 50))["strength_loy"]

## PP
fit <- feols(pp_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
             data = subset(data_final, eus==1), vcov=conley_vcov)
conleys[3,1] <- coef(fit)["strength_loy"]
# 30 km (base)
conleys[3,2] <- se(fit)["strength_loy"]
# 40 km
conleys[3,3] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 40))["strength_loy"]
# 50 km
conleys[3,4] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 50))["strength_loy"]

## PSOE
fit <- feols(psoe_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
             data = subset(data_final, eus==1), vcov=conley_vcov)
conleys[4,1] <- coef(fit)["strength_loy"]
# 30 km (base)
conleys[4,2] <- se(fit)["strength_loy"]
# 40 km
conleys[4,3] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 40))["strength_loy"]
# 50 km
conleys[4,4] <- se(fit, vcov=vcov_conley(lat="lat_rec", lon = "lon_rec", cutoff = 50))["strength_loy"]

colnames(conleys) <- c("Coef.", "S.E. (30 km)", "S.E. (40 km)", "S.E. (50 km)")
rownames(conleys) <- c("Rad. Nat. (Avg)", "PNV (Avg)", "PP (Avg)", "PSOE (Avg)")

## Export table
stargazer(conleys, out = paste0(outtabledir, "tableC5.tex"))


##### Table C6: WCBSE #####
## Matrix to store the results
rob <- matrix(nrow=4, ncol=2)

fit <- feols(ia_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
             data = subset(data_final, eus==1))
rob[1,1] <- coef(fit)["strength_loy"]; rob[1,2] <- boottest(fit, B=999, param="strength_loy", clustid = "com", seed = 123)$p_val

fit <- update(fit, pnv_share_avg ~ .)
rob[2,1] <- coef(fit)["strength_loy"]; rob[2,2] <- boottest(fit, B=999, param="strength_loy", clustid = "com", seed = 123)$p_val

fit <- update(fit, pp_share_avg ~ .)
rob[3,1] <- coef(fit)["strength_loy"]; rob[3,2] <- boottest(fit, B=999, param="strength_loy", clustid = "com", seed = 123)$p_val

fit <- update(fit, psoe_share_avg ~ .)
rob[4,1] <- coef(fit)["strength_loy"]; rob[4,2] <- boottest(fit, B=999, param="strength_loy", clustid = "com", seed = 123)$p_val

colnames(rob) <- c("Coef.", "p-value")
rownames(rob) <- c("Rad. Nat. (Avg)", "PNV (Avg)", "PP (Avg)", "PSOE (Avg)")


## Export table
stargazer(rob, out = paste0(outtabledir, "tableC6.tex"))


##### Figure C1: Trimming extreme values #####

# Empty matrix to store results
range <- c(0.01, 0.02, 0.03, 0.04, 0.05)

p1 <- ggplot(do_trim("ia_share_avg"), aes(x=coef, y=as.factor(qt))) + geom_point() + geom_errorbar(aes(xmin=lower, xmax=upper)) +
  geom_vline(xintercept=0, linetype="dashed", color="blue") + labs(x="Coeff. signal strength", y="Prop. trimmed", title="Radical Nationalists") +
  theme_bw(base_size=9)

p2 <- ggplot(do_trim("pnv_share_avg"), aes(x=coef, y=as.factor(qt))) + geom_point() + geom_errorbar(aes(xmin=lower, xmax=upper)) +
  geom_vline(xintercept=0, linetype="dashed", color="blue") + labs(x="Coeff. signal strength", y="Prop. trimmed", title="PNV") +
  theme_bw(base_size=9)

p3 <- ggplot(do_trim("pp_share_avg"), aes(x=coef, y=as.factor(qt))) + geom_point() + geom_errorbar(aes(xmin=lower, xmax=upper)) +
  geom_vline(xintercept=0, linetype="dashed", color="blue") + labs(x="Coeff. signal strength", y="Prop. trimmed", title="PP") +
  theme_bw(base_size=9)

p4 <- ggplot(do_trim("psoe_share_avg"), aes(x=coef, y=as.factor(qt))) + geom_point() + geom_errorbar(aes(xmin=lower, xmax=upper)) +
  geom_vline(xintercept=0, linetype="dashed", color="blue") + labs(x="Coeff. signal strength", y="Prop. trimmed", title="PSOE") +
  theme_bw(base_size=9)

ggarrange(p1, p2, p3, p4, nrow = 2, ncol=2)
ggsave(filename = "figureC1.pdf",
       path = outplotdir, height = 8, units = "cm")


##### Figure C2: Leave-one-out tests for influential observations #####

## IA 

### Loop for leave-one-out ###
loo_loy <- matrix(, nrow = nrow(data_final[data_final$eus==1,]), ncol = 4)


for (i in 1:nrow(data_final[data_final$eus==1,])) {
  # Keep only the BC sample
  d <- data_final %>% filter(eus==1)
  # Drop town i
  d <- d[-i,]
  
  # Estimate the model for IA in the data without town i
  fit <- feols(ia_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
               data = d, vcov = conley_vcov)
  
  # Extract the coefficient
  coef <- coef(fit)["strength_loy"]
  
  # Compute CI
  ci <- confint(fit)["strength_loy",]
  
  # Store the coefficients and their CI in the matrix
  loo_loy[i, 1] <- i
  loo_loy[i, 2] <- coef
  loo_loy[i, 3] <- as.numeric(ci[1])
  loo_loy[i, 4] <- as.numeric(ci[2])
  
  
}



## Plot of LOO coefficients and CIs
pl1 <- loo_loy %>% 
  as.data.frame() %>%
  ggplot(aes(x=V1, y=V2)) +
  geom_point(color = "red") + geom_errorbar(aes(ymin = V3, ymax = V4), color="grey") +
  geom_hline(yintercept = 0, color="blue", linetype = "dashed") +
  scale_x_discrete(labels = NULL) +
  labs(x = "", y = "Coeff. signal strength", title="Radical Nationalists") +
  theme_classic()



## PNV
### Loop for leave-one-out ###
loo_loy <- matrix(, nrow = nrow(data_final[data_final$eus==1,]), ncol = 4)


for (i in 1:nrow(data_final[data_final$eus==1,])) {
  # Keep only the BC sample
  d <- data_final %>% filter(eus==1)
  # Drop town i
  d <- d[-i,]
  
  # Estimate the model for PNV in the data without town i
  fit <- feols(pnv_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
               data = d, vcov = conley_vcov)
  
  # Extract the coefficient
  coef <- coef(fit)["strength_loy"]
  
  # Compute CI
  ci <- confint(fit)["strength_loy",]
  
  # Store the coefficients and their CI in the matrix
  loo_loy[i, 1] <- i
  loo_loy[i, 2] <- coef
  loo_loy[i, 3] <- as.numeric(ci[1])
  loo_loy[i, 4] <- as.numeric(ci[2])
  
  
}



## Plot of LOO coefficients and CIs
pl2 <- loo_loy %>% 
  as.data.frame() %>%
  ggplot(aes(x=V1, y=V2)) +
  geom_point(color = "red") + geom_errorbar(aes(ymin = V3, ymax = V4), color="grey") +
  geom_hline(yintercept = 0, color="blue", linetype = "dashed") +
  scale_x_discrete(labels = NULL) +
  labs(x = "", y = "Coeff. signal strength", title="PNV") +
  theme_classic()


## PP
### Loop for leave-one-out ###
loo_loy <- matrix(, nrow = nrow(data_final[data_final$eus==1,]), ncol = 4)


for (i in 1:nrow(data_final[data_final$eus==1,])) {
  # Keep only the BC sample
  d <- data_final %>% filter(eus==1)
  # Drop town i
  d <- d[-i,]
  
  # Estimate the model for PP in the data without town i
  fit <- feols(pp_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
               data = d, vcov = conley_vcov)
  
  # Extract the coefficient
  coef <- coef(fit)["strength_loy"]
  
  # Compute CI
  ci <- confint(fit)["strength_loy",]
  
  # Store the coefficients and their CI in the matrix
  loo_loy[i, 1] <- i
  loo_loy[i, 2] <- coef
  loo_loy[i, 3] <- as.numeric(ci[1])
  loo_loy[i, 4] <- as.numeric(ci[2])
  
  
}



## Plot of LOO coefficients and CIs
pl3 <- loo_loy %>% 
  as.data.frame() %>%
  ggplot(aes(x=V1, y=V2)) +
  geom_point(color = "red") + geom_errorbar(aes(ymin = V3, ymax = V4), color="grey") +
  geom_hline(yintercept = 0, color="blue", linetype = "dashed") +
  scale_x_discrete(labels = NULL) +
  labs(x = "", y = "Coeff. signal strength", title="PP") +
  theme_classic()



## PSOE
### Loop for leave-one-out ###
loo_loy <- matrix(, nrow = nrow(data_final[data_final$eus==1,]), ncol = 4)


for (i in 1:nrow(data_final[data_final$eus==1,])) {
  # Keep only the BC sample
  d <- data_final %>% filter(eus==1)
  # Drop town i
  d <- d[-i,]
  
  # Estimate the model for PSOE in the data without town i
  fit <- feols(psoe_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
               data = d, vcov = conley_vcov)
  
  # Extract the coefficient
  coef <- coef(fit)["strength_loy"]
  
  # Compute CI
  ci <- confint(fit)["strength_loy",]
  
  # Store the coefficients and their CI in the matrix
  loo_loy[i, 1] <- i
  loo_loy[i, 2] <- coef
  loo_loy[i, 3] <- as.numeric(ci[1])
  loo_loy[i, 4] <- as.numeric(ci[2])
  
  
}



## Plot of LOO coefficients and CIs
pl4 <- loo_loy %>% 
  as.data.frame() %>%
  ggplot(aes(x=V1, y=V2)) +
  geom_point(color = "red") + geom_errorbar(aes(ymin = V3, ymax = V4), color="grey") +
  geom_hline(yintercept = 0, color="blue", linetype = "dashed") +
  scale_x_discrete(labels = NULL) +
  labs(x = "", y = "Coeff. signal strength", title="PSOE") +
  theme_classic()


ggarrange(pl1, pl2, pl3, pl4, nrow=4, ncol=1)

## Export the plot
ggsave(filename = "figureC2.png", width=12, height=16, units="cm",
       path = outplotdir)


################################################################################
##### Interactive models #####
################################################################################

##### Table 4: Pre-war nationalism #####
int <- NULL

## IA
int[[1]] <- feols(ia_share_avg ~ strength_loy*pnv_share_33 + fs_loy + log_height + log_area + log_pop_70 | prov,
                  data = subset(data_final, eus==1),
                  vcov = conley_vcov)

## PNV
int[[2]] <- update(int[[1]], pnv_share_avg ~ .)

## PP
int[[3]] <- update(int[[1]], pp_share_avg ~ .)

## PSOE
int[[4]] <- update(int[[1]], psoe_share_avg ~ .)

## Export table
etable(int, keep = c("Signal strength", "PNV"), tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "table4.tex"), replace = TRUE)



##### Table 5: Share of Basque speakers #####
int <- NULL

## IA
int[[1]] <- feols(ia_share_avg ~ strength_loy*euskaldun_share + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                  data = subset(data_final, eus==1),
                  vcov = conley_vcov)

## PNV
int[[2]] <- update(int[[1]], pnv_share_avg ~ .)

## PP
int[[3]] <- update(int[[1]], pp_share_avg ~ .)

## PSOE
int[[4]] <- update(int[[1]], psoe_share_avg ~ .)

## Export table
etable(int, keep = c("Signal strength", "Basque"), tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "table5.tex"), replace = TRUE)



##### Figure C3: Binning estimators with PNV v.s. #####
## IA
set.seed(123)
if1 <- interflex(Y = "ia_share_avg", D = "strength_loy", X = "pnv_share_33", 
                 Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33"), FE = "prov", 
                 data = subset(data_final, eus==1 & !is.na(ia_share_avg)), 
                 estimator = "binning", cl = "com", vartype = "bootstrap", 
                 bin.labs = F, ylab = "Marginal effect", xlab = "Moderator: PNV '33", subtitles = "Rad. Nat. (Avg.)", theme.bw = TRUE,
                 parallel=F)

## PNV
set.seed(123)
if2 <- interflex(Y = "pnv_share_avg", D = "strength_loy", X = "pnv_share_33", 
                 Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33"), FE = "prov", 
                 data = subset(data_final, eus==1 & !is.na(pnv_share_avg)), 
                 estimator = "binning", cl = "com", vartype = "bootstrap", 
                 bin.labs = F, ylab = "Marginal effect", xlab = "Moderator: PNV '33", subtitles = "PNV (Avg.)", theme.bw = TRUE,
                 parallel=F)

## PP
set.seed(123)
if3 <- interflex(Y = "pp_share_avg", D = "strength_loy", X = "pnv_share_33", 
                 Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33"), FE = "prov", 
                 data = subset(data_final, eus==1 & !is.na(pp_share_avg)), 
                 estimator = "binning", cl = "com", vartype = "bootstrap", 
                 bin.labs = F, ylab = "Marginal effect", xlab = "Moderator: PNV '33", subtitles = "PP (Avg.)", theme.bw = TRUE,
                 parallel=F)

## PSOE
set.seed(123)
if4 <- interflex(Y = "psoe_share_avg", D = "strength_loy", X = "pnv_share_33", 
                 Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33"), FE = "prov", 
                 data = subset(data_final, eus==1 & !is.na(psoe_share_avg)), 
                 estimator = "binning", cl = "com", vartype = "bootstrap", 
                 bin.labs = F, ylab = "Marginal effect", xlab = "Moderator: PNV '33", subtitles = "PSOE (Avg.)", theme.bw = TRUE,
                 parallel=F)

## Plot
ggarrange(if1$figure, if2$figure, if3$figure, if4$figure, ncol=2, nrow = 2)


## Export the plot
ggsave(filename = "figureC3.pdf", path = outplotdir,width = 14,  height = 14, units = "cm")



################################################################################
##### Table 3: Other measures of mobilization #####
################################################################################

models <- NULL

## Turnout referendum 1978
models[[1]] <- feols(turnout_78 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                     data = subset(data_final, eus==1),
                     vcov = conley_vcov)

## No share referendum 1978
models[[2]] <- update(models[[1]], no_share_78 ~ .)

## Turnout referendum 1979
models[[3]] <- update(models[[1]], turnout_r79 ~ .)

## Yes share referendum 1979
models[[4]] <- update(models[[1]], yes_share_r79 ~ .)


## Export table
etable(models, keep = c("Signal strength"), tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "table3.tex"), replace = TRUE)





################################################################################
##### Language adoption #####
################################################################################

##### Figure G1: Basque speakers in 1972 and change over time #####
data_final %>% filter(eus==1) %>% 
  ggplot(aes(x=euskaldun_share, y=change_72_16)) +
  geom_point() + geom_smooth(method = "loess") + 
  labs(x="Euskera share 1972", y = "Euskera share diff. (1972-2016)") + theme_bw()

ggsave(filename = "figureG1.pdf", height = 10, unit = "cm", path = outplotdir)


##### Figure 5: Regression model through interflex #####
set.seed(1234)
interflex(Y="change_72_16", D="strength_loy", X="euskaldun_share",
          Z=c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33", "euskaldun_share"), FE = "prov",
          data = subset(data_final, eus==1 & !is.na(change_72_16)),
          estimator = "binning", cl="com", vartype = "bootstrap",
          bin.labs = F, ylab = "Marginal effect", xlab = "Moderator: Euskera share 1972", theme.bw = T,
          parallel=F)

ggsave(filename = "figure5.pdf", height = 10, unit = "cm", path = outplotdir)



################################################################################
##### Figure D1: Radio and radical nationalism in Navarra #####
################################################################################

## Note in this sample there is only one province, so no province FE
nav <- NULL
nav[[1]] <- feols(IA_share_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33,
                  data = subset(data_final, eus==0),
                  vcov = conley_vcov)
nav[[2]] <- update(nav[[1]], hb_share_79 ~ .)
nav[[3]] <- update(nav[[1]], hb_share_82 ~ .)
nav[[4]] <- update(nav[[1]], hb_share_86 ~ .)
nav[[5]] <- update(nav[[1]], hb_share_89 ~ .)
nav[[6]] <- update(nav[[1]], hb_share_93 ~ .)
nav[[7]] <- update(nav[[1]], hb_share_96 ~ .)
nav[[8]] <- update(nav[[1]], amaiur_share_11 ~ .)
nav[[9]] <- update(nav[[1]], bildu_share_15 ~ .)
nav[[10]] <- update(nav[[1]], bildu_share_16 ~ .)
nav[[11]] <- update(nav[[1]], bildu_share_19 ~ .)
nav[[12]] <- update(nav[[1]], ia_share_avg ~ .)

## Plot and save
cbind(c(1:length(nav)),
      as.factor(c(election_years[!election_years%in%c("2000", "2004", "2008")], "Average")),
      do.call("rbind", lapply(nav, function(x) coef(x)["strength_loy"])),
      do.call("rbind", lapply(nav, function(x) confint(x)["strength_loy",]))) %>%
  as.data.frame() %>%
  dplyr::rename(year = `c(1:length(nav))`, coefs = strength_loy, lower = `2.5 %`, upper = `97.5 %`) %>%
  ggplot(aes(x=year, y=coefs)) + geom_point() + geom_errorbar(aes(ymin=lower, ymax=upper), color="black") +
  geom_hline(yintercept=0, linetype="dashed", color="blue") +
  scale_x_discrete(limits = factor(1:length(nav)), labels = c(election_years[!election_years%in%c("2000", "2004", "2008")], "Average"), drop=F) +
  labs(x="", y="", title="Radical nationalists, Navarre") + theme_bw() + theme(axis.text.x = element_text(size=8))

ggsave(filename = "figureD1.pdf", path = outplotdir, width = 12, units = "cm")



################################################################################
##### Table C7: Results with regional elections #####
################################################################################

### Outcomes (BC) ###
## PNV: 1980-2020
## IA: HB (1980-1994), EH (1998-2001), EH Bildu (2012-2020)
## PP: AP/CP (1980-1986), PP (1990-2020)
## PSOE: 1980-2020

## Regressions on average outcomes
reg_avgs <- NULL

## IA
reg_avgs[[1]] <- feols(ia_er_share_avg ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                       data=subset(data_final, eus==1),
                       vcov=conley_vcov)

## PNV
reg_avgs[[2]] <- update(reg_avgs[[1]], pnv_er_share_avg ~ .)

##PP
reg_avgs[[3]] <- update(reg_avgs[[1]], pp_er_share_avg ~ .)

## PSOE
reg_avgs[[4]] <- update(reg_avgs[[1]], psoe_er_share_avg ~ .)

## Export table
etable(reg_avgs, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC7.tex"), replace = TRUE)


################################################################################
##### Table C8: Results on UCD (1979) #####
################################################################################
mod_ucd <- feols(ucd_share_79 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                 data = subset(data_final, eus==1),
                 vcov = conley_vcov)

## Export table
etable(mod_ucd, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC8.tex"), replace = TRUE)


################################################################################
##### Combined radio signals #####
################################################################################


##### Table C9: Average of RPL and RPS #####
lt_avgs_gip <- NULL

## IA
lt_avgs_gip[[1]] <- feols(ia_share_avg ~ avg_strength_gip + avg_freespacestrength_gip + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                           data = subset(data_final, eus==1),
                           vcov = conley_vcov)

## PNV
lt_avgs_gip[[2]] <- update(lt_avgs_gip[[1]], pnv_share_avg ~ .)

## PP
lt_avgs_gip[[3]] <- update(lt_avgs_gip[[1]], pp_share_avg ~ .)

## PSOE
lt_avgs_gip[[4]] <- update(lt_avgs_gip[[1]], psoe_share_avg ~ .)

## Export table
etable(lt_avgs_gip, keep = "Avg", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC9.tex"), replace = TRUE)


##### Table C10: RPL and RPS enter separately #####
lt_avgs_sep <- NULL

## IA
lt_avgs_sep[[1]] <- feols(ia_share_avg ~ strength_loy + strength_rpss + fs_loy + freespacestrength_rpss + log_height + log_area + log_pop_70 + pnv_share_33 | prov,
                          data = subset(data_final, eus==1),
                          vcov = conley_vcov)

## PNV
lt_avgs_sep[[2]] <- update(lt_avgs_sep[[1]], pnv_share_avg ~ .)

## PP
lt_avgs_sep[[3]] <- update(lt_avgs_sep[[1]], pp_share_avg ~ .)

## PSOE
lt_avgs_sep[[4]] <- update(lt_avgs_sep[[1]], psoe_share_avg ~ .)

## Export table
etable(lt_avgs_sep, keep = "Signal strength", tex= T, style.tex = style.tex("aer"), file = paste0(outtabledir, "tableC10.tex"), replace = TRUE)


################################################################################
##### Figure 6: ETA activities #####
################################################################################


##### ETA attacks #####
attacks <- NULL
## 1968-1977
attacks[[1]] <- feols(log_eta_att_68_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
attacks[[2]] <- feols(any_eta_att_68_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
## 1978-1994
attacks[[3]] <- feols(log_eta_att_78_94 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
attacks[[4]] <- feols(any_eta_att_78_94 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
## 1995-2003
attacks[[5]] <- feols(log_eta_att_95_03 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
attacks[[6]] <- feols(any_eta_att_95_03 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)


##### ETA victims #####
victims <- NULL
## 1968-1977
victims[[1]] <- feols(log_eta_vict_68_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
victims[[2]] <- feols(any_eta_vict_68_77 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
## 1978-1994
victims[[3]] <- feols(log_eta_vict_78_94 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
victims[[4]] <- feols(any_eta_vict_78_94 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
## 1995-2003
victims[[5]] <- feols(log_eta_vict_95_03 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
victims[[6]] <- feols(any_eta_vict_95_03 ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)

##### ETA prisoners #####
prisoners <- NULL
prisoners[[1]] <- feols(log_pris ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)
prisoners[[2]] <- feols(any_pris ~ strength_loy + fs_loy + log_height + log_area + log_pop_70 + pnv_share_33 | prov, data = subset(data_final, eus==1), vcov = conley_vcov)

##### Plot #####
ggarrange(plot_eta_coefs(attacks) + labs(subtitle="Attacks"),
          plot_eta_coefs(victims) + labs(subtitle="Victims") + theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), axis.title.y=element_text(angle=0)),
          plot_eta_coefs2(prisoners) + labs(subtitle="Prisoners") + theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), axis.title.y=element_text(angle=0)),
          ncol=2, nrow=2, common.legend=T, legend="right", align="hv")

ggsave(filename=paste0(outplotdir, "figure6.pdf"), width=16, height=10, units="cm")



################################################################################
##### Figure E1: Associations #####
################################################################################
## Generate the association density per 100 inhab. (using 1970 population)
data_final$assoc_tot_pc <- (data_final$assoc_tot/data_final$pop_70)*100
data_final$assoc_cult_pc <- (data_final$assoc_cult/data_final$pop_70)*100
data_final$assoc_cult2_pc <- (data_final$assoc_cult2/data_final$pop_70)*100

set.seed(1234)
assoc1 <- interflex(Y = "assoc_tot_pc", D = "strength_loy", X = "euskaldun_share", 
          Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33", "euskaldun_share"), FE = "prov", 
          data = subset(data_final, eus==1 & !is.na(euskaldun_share) & !is.na(assoc_tot_pc)), 
          estimator = "binning", cl = "com", vartype = "bootstrap", 
          bin.labs = F, ylab = "", xlab = "", subtitles = "All associations", theme.bw = TRUE,
          parallel=F)

set.seed(1234)
assoc2 <- interflex(Y = "assoc_cult_pc", D = "strength_loy", X = "euskaldun_share", 
          Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33", "euskaldun_share"), FE = "prov", 
          data = subset(data_final, eus==1 & !is.na(assoc_cult_pc) & !is.na(euskaldun_share)), 
          estimator = "binning", cl = "com", vartype = "bootstrap", 
          bin.labs = F, ylab = "", xlab = "", subtitles = "Cultural (broad)", theme.bw = TRUE,
          parallel=F)

set.seed(1234)
assoc3 <- interflex(Y = "assoc_cult2_pc", D = "strength_loy", X = "euskaldun_share", 
          Z = c("fs_loy", "log_height", "log_area", "log_pop_70", "pnv_share_33", "euskaldun_share"), FE = "prov", 
          data = subset(data_final, eus==1 & !is.na(assoc_cult2_pc) & !is.na(euskaldun_share)), 
          estimator = "binning", cl = "com", vartype = "bootstrap", 
          bin.labs = F, ylab = "", xlab = "", subtitles = "Cultural (specific)", theme.bw = TRUE,
          parallel=F)

## Plot
ggarrange(assoc1$figure, assoc2$figure, assoc3$figure, ncol=3, nrow = 1) %>%
  annotate_figure(left = "Marginal effect", bottom = "Moderator: Euskera share 1972")


## Export the plot
ggsave(filename = "figureE1.pdf", path = outplotdir,width = 24, height=10, units = "cm")





